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Summary. We study the optimal design of switching measurements of small Josephson junc- 
tion circuits which operate in the macroscopic quantum tunnelling regime. Starting from the 
D-optimality criterion we derive the optimal design for the estimation of the unknown parame- 
ters of the underlying Gumbel type distribution. As a practical method for the measurements, 
we propose a sequential design that combines heuristic search for initial estimates and max- 
imum likelihood estimation. The presented design has immediate applications in the area of 
superconducting electronics implying faster data acquisition. The presented experimental re- 
sults confirm the usefulness of the method. 
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1. Introduction 

Randomness plays a central role in quantum physics. On the fundamental level, randomness 
is not caused only by unknown cova riates or measure ment error but is understood to be 
an intrinsic property of Nature itself <)Ballentinelll998(l . From the statistical point of view, 
this leads to fascinating problems where the distributional assumptions arise directly from 
the fundamental laws of physics. 

The Josephson junctions (JJ) are important non-linear components of superconduct- 
ing electronics. A typical junction involves two superconducting electrodes separated by 
an insulator gap l)JosephsorI I19 62). In practice, miniature J J circuits are manufactured 
by lithographic means and operated at cryogenic temperatures, i.e., below 4.2 K. The 
strong dependence of the physical parameters of JJ circuits as function of changes in en- 
vironmental variables, for instance, temperature, electric noise, and magnetic field make s 
the JJ circuits to have several applications as ultra-sensitive sensors l|Pekola et all 12005). 
Moreover, certain JJ circuits a re promising candidates for realization of quantum com- 
putation dMakhlin et all l200ll). F or a thorough review of JJ circuits see, for instance, 
l|Grabert and Devoret fEds.llll99'^ . 

An experiment called switching measurement, see Fig.^ is a common way to probe the 
properties of a JJ circuit sample. In the experiment, sequences of current pulses are applied 
to the sample, while the voltage over the structure is monitored. The superconducting 
electrodes of the sample are described by a complex valued quantum mechanical wave 
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function and the physical state of the JJ is described by the phase difference over the 
junction. The dynamics of the phase difference arc analogous to those of a massive particle 
in a tilted cosine potential, see Fig. ^d). A sufficiently high current pulse causes the phase 
difference to escape from the local minimum of the potential to the direction of the steepest 
descent. In the experiment the observed voltage pulse indicates that switching has occurred. 
The appearance of a voltage pulse to a single applied current pulse, being governed by the 
laws of quantum mechanics, is purely random. Below, we call observed voltage pulse by 
Y = 1, otherwise Y = 0. The analysis of the statistics of this random variable with respect 
to current pulse heights yields information about the physical parameters of a J J circuit. 

We study the switching dynamics of a single J J. Quantum mechanical arguments studied 
in more detail in Section |21 suggest that the switching rate is of the form 7 = e ax+b , where 
x = x (I) is a monotonic function depending on the height / of the applied current pulse and 
a and b are the parameters depending on the geometry and environment of the junction. 
Moreover, we note that x{I) must be independent of parameters a and b. Assuming constant 
switching rate during the pulse we get the probability of measuring a voltage pulse at an 
applied current pulse 

P(Y = 1) = 1 — e - cx p( ax + b )_ (1) 

For the simplicity we consider the case x(I) = I in this paper. In statistics, model J3> is 
often called complementary log-log regression. 

Assume that a sequence of n current pulses of height i — 1, 2, . . . , n is applied to 
the sample. The obtained data consist of covariate values Xi, i — 1,2, ...,n and cor- 
responding binary responses yt. Our aim is to choose the covariates such that the pa- 
rameters a and b can be estimated as well as p ossible from the measured data. In sta- 
tistical terms this is an optimal design problem ijAtkinson and DonevL Il992t IPukelsheiml 
Il993l iLiski et all l2002t iBereer and WoneL |2005|) . The literature on optimal designs is 
concentrated on linear models but work has been done also with nonlinear models . Sev- 
eral authors fealish and Rosenbergerl Il978t lAbdelbasit and PlackettL Il98i iMinkinl. Il987t 
ISitter and WuUl993llMathew and SinhaLl200l|) have considered optimal designs for binary 
data in the case of logistic regression 

p ( Y = D = 1 + e -U - (2) 

where a and b are unknown parameters and x is a covariate. Compared to the logistic 
regression the complementary log-log regression Q is less frequently studied or applied. 
Complementary log-log regression is often used as an alternative model for logistic regres- 
sion @ if the symmetry assumption does not fit to the data, but in our application the 
complementary log-log regression arises directly from the physical properties of JJs. The 
ke y reference c onsid ering optimal design for complementary log-log regression is the work 
bv lFord et al . (1992). They provide the general framework for the D and c optimal designs 
and the solution for complementary log-log regression is obtained as a special case of their 
general results. 

In this paper we propose a sequential two-point design for estimating the unknown pa- 
rameters of a JJ circuit. The design is derived using D-optimality criterion, which equals to 
minimizing the generalized variance of the parameters. Based on the data available, the de- 
sign yields the optimal current pulse heights for the next measurement in real time. We are 
ultimately interested in minimizing the time needed for an experiment. The proposed se- 
quential design also takes into account the special requirements of switching measurements, 
e.g. the time needed for the reconfiguration of the signal generator. 
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Fig. 1. Illustration of a switching measurement, (a) Scanning electron micrograph of a typical two- 
junction JJ sample, (b) Sketch of applied current (I) pulses and resulting (V) voltage pulses, (c) 
The associated switching probability as a function of the height of the current pulses, (d) Quantum 
mechanical process leading to switching: The phase "particle" is initially localized in a local minimum 
of the potential. When the current pulse is applied, the phase particle tunnels through the potential 
barrier and escapes downhill. 
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The rest of the paper is organized as follows: In Section [3 a more detailed description 
of physical background of the switching measurements is given. In Section |3J the optimal 
design for the complementary log-log regression defined in equation JQl is derived. In 
Section 21 we propose a sequential design for switching measurements. In Section [5] the 
sequential design is applied to switching measurements of a J J circuit. Results from the 
real data are compared with simulation results. Finally, Section concludes the paper. 

2. Physical background 

In this section we study the quantum mechanical arguments leading to the model . The 
switching dynamics are described by a Poisson process where the switching rate 7 is constant 
over the current pulse which has duration At. The probability of switching is then obtained 
as 

P(Y = 1) = 1 - e- Atl . (3) 

The duration At of the current pulse is very short and remains fixed in the experiment. 
Thus, we consider At as a part of the switching rate and do not write it explicitly in the 
subsequent equations. The relation between the switching rate 7 and the height of the 
current pulse / can be written in the form 

7 (7) = e ax{I)+ \ (4) 

where x{I) is a monotonic function of the bias current / and a and b are constants depending 
on the physical parameters of the J J circuit. A commonly used approximation is to assume 
that the activation threshold amplitude is linearly proportional to current I. This leads to 
the switching rate 

1(1) = e aI+b . (5) 

Although not based on microscopic theory, the switching rate model (0 describes the results 
of typical switching measurements quite well, see Fig. 

In order to obtain the functional form of x(I) in the switching rate model we need 
a more accurate physical model to describe the system under consideration. Theoretically, 
the superconducting electrodes of the JJ sample are modeled by a complex valued quantum 
mechanical wave function and the physical state of the J J is described by the phase difference 
over the junction. The dynamics of the phase difference (f> are analogous to those of a massive 
particle in a tilted cosine potential U under gravitational force, 

U = -Ej(— cf> + cos cp), (6) 

where Ej is so-called Josephson energy, / is the magnitude of the applied bias current, and 
Ic is the critical current for which the potential changes into monotonic function. The 
potential is depicted during a pulse in Fig.QJd). The applied current changes the shape of 
the potential landscape of a JJ circuit in a controlled way. The phase difference <f> initially 
set into one of the local minima becomes a metastable state for finite current /; The state 
may decay by escape of the phase from the local minimum of the potential to the direction 
of the steepest descent. This escape event is called switching. In the experiment the voltage 
over the JJ sample is monitored and the observed voltage pulse indicates that switching has 
occurred. 
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Fig. 2. Switching rate 7 as a function of the relative bias current I/Ic- The solid line represents the 
switching rate 7mqt in the macroscopic quantum tunnelling, while the dotted and the dashed line rep- 
resents the switching rate 7ta under the thermal activation model in the temperatures 1 K and 0.1 K, 
respectively. The JJ parameters are set to have typical values: Josephson energy Ej = 1CT 21 J and 
capacitance C = 50 femtofarads. 7mqt is calculated on the high Q limit. The relation between the 
switching rate 7 and the height of the current pulse / is nearly linear on the logarithmic scale, which 
supports the use of the approximation -y(I) = e aI+b . 
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Depending on the temperature, the model ^ is explained by the thermal activation 
model of classical physics or by quantum mechanical analysis. In an experiment in tem- 
perature above 0.2 K, there will be thermal fluctuations in the bias current / as the JJ 
is connected to the measurement wiring. These fluctuations make the phase particle to 
absorb energy and eventually escape from the potential well. One of the central results in 
statistical physics is that the probability for system to be in energy state E at temperature 
T is proportional to exp(— E/U-qT), where fee ~ 1.381 • 10~ 23 J/K is the Bolztman constant. 
In order to escape from the well the phase particle needs to have more energy than the 
barrier height AC/. Based on these principles the average escape rate for the potential of 
equation © can be shown to follow the Kramer's formula 

7ta(/) = exp(-Af/(/)/fc B T), (7) 

where AU is the height of the potential barrier 

w m = ^ (l -A.Y<\ (8) 



3 V Ic 
and cup is the characteristic oscillation frequency 



E 3 e 



Ic 



where C is the capacitance of the JJ, e = 1.602 ■ 10 -19 coulombs is the elementary charge 
and h « 1.055 • 10~ 34 Js is the reduced Planck's constant. In above equations Ej and C 
are the unknown quantities - Ic is a function of Ej. It is easy to see that the switching 
rate J7|) can be written in the form 7ta(-0 = e ax ^ +b . 

When a JJ is cooled to tens of millike lvins one can still observe the phase part icle es- 
caping from the quasistationary potential l|Voss and W ebb. 1981t lJackel et allll98lf) . Since 
the thermal activation, equation Q yields vanishing switching rate for that temperature 
range, the microscopic mechanism must be differe nt from the one described a bove. This 
process is called macroscopic quantum tunnelling IjCaldeira and Leggettl Il980j) and to ex- 
plain it^d^taikedquantum mechanical analysis is required. The analysis leads to switching 
rate l|Grabert et all Il985|) 

1mqt{I) = A(I)e- B ^ . (10) 

The prefactor in equation l(TT)|) is 



where x( a ) = 12v / 67r(l + ca + 0(a 2 ), with c « 2.86, and Q is the quality factor of the 
oscillator describing the effective environment of the junction. The exponent reads 



where s(a) = 36/5[l+45C(3)/7r 3 a+C(a 2 )] and ((3) » 1.202 is the value of the Riemann zeta 
function at 3. Again, the switching rate (fTT)f) can be written in the form 7mqt(-0 = 



e ax(I)+b 
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3. D-optimal designs for complementary log-log regression 



Assume that a binary random variable Y is distributed as follows 



P(Y = 1) 



1 



- cxp(aa;-f b) 



_ n \ _ - exp(ax+b) 



P{Y = 0) = e 



(13) 



where a and b are unknown parameters and x is a covariate. The cumulative distribution 
function 



F(z) = 1 



- cxp(2:) 



(14) 



where z — ax + b, is known in literature by names Gumbel distribution, Gompertz distribu- 
tion, LogWeibull distribution, Fisher-Tippett distribution and extreme value distribution. 
Whereas, the inverse of (|14H 



F 1 (u) = ax + b = log(- log(l - u)) 



(15) 



is called complementary log-log link and is one of the standard link functions in generalized 
linear models. 

The problem is to select such covariate values x±, %2, . . . ,x n that the parameters a and 
b can be estimated from the observed data (xi,yi), i = 1, 2, . . . , n in optimal manner. Here 
we define the optimality as D-optimality, which equals to maximizing the determinant of the 
information matrix, and estimate the parameters by maximum likelihood. Starting from 
the definition (fHfl) we obtain the log-likelihood of the data (xi,yi), i = 1, 2, . . . , n 



L(a, b) = J2 (Vi log(l - e" - (1 - Vi )e 

The information matrix 



axi~\-b 



takes the form 



dbda db 2 

, n n v 



(16) 



(17) 



n 



where = ax, + b and 



i=l 

2z 



e exp(z) _ ' 

The D-optimality is equivalent to finding such xi,x%, ■ ■ • ,x n that 



(18) 



(19) 



det(J) 



i=l J Li=l 

n n — 1 

^29(zi)g(zj)(xi-x 

j— z+1 1 



(20) 
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is m aximized . It is assu med in the ma ximization that the true values of a and b are known 
(see l|Minkinl.ll987|) and l|Sitterl Il99^ for a discussion). 

First we will find a locally optimal two-point design n = 2. This simplifies equation 120|) 

to 

det(J) = \g(z 1 )g(z 2 )(z 1 - z 2 f . (21) 

The contour plot of det(J) as a function z\ and z 2 is depicted in Fig.|3 The values z\ and 
z 2 that maximize det(J) 

z\ « 0.979633 
z* 2 w -1.33774 

(22) 

can be solved numerically from l|21|l . The corresponding response probabilities are 

F[z\) w 0.930295 

F{z* 2 ) w 0.230826, (23) 

indicating that the design is asymmetric with respect to the middle point of the response 
curve. Since the function is symmetric with respect to exchanging of the variables, the 
mirror points provide also the same maximum value. 

In order to show that the two-point maximum (122(1 is a local maximum also when n > 2, 
we define functions 

v(zi, zi) = g(z 1 )g(z 2 )(z 1 - z 2 ) 2 (24) 
d 

V (Zi,Z2) = —v(zi,Z2). (25) 

By inspecting the functional form of v(z{, Zj) we observe that 

v(z*, z|) = v{z~2, Zi) => local maximum of v (26) 

v{z, z) = local minimum of v. (27) 

We are searching for a local maximum of 

-j n n—1 

det(J) = - Y, E^'^)- ( 28 ) 



a 

j—i+l 2—1 



For a local extreme of dct(J) it holds 



v(z h Zj)=0 for alii. (29) 



Because v (z*, z|) = 0, a solution for equation Ij29(l is obtained if 

Zi£{zl,z%} foralH (30) 
From equation Q26J) we conclude that a local maximum is obtained, for instance, when 

Zi = z\ if i is odd 

Zi — z 2 Hi is even. (31) 
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Fig. 3. Contours of the ratio 100 ■ det(J(2i, 2 2 ))/det(J(«i, zl)), where det(J(zi, z 2 )) is the determi- 
nant of the information matrix for the two-point design [z\,zi). The maximum of det(J) is achieved 
at the point (z* « 0.98, zl « -1.34) and at its mirror point. 
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For the co mplete proof that the solution IL'21) is th e global maximum, we r efer to the general 
result s bv lFord et all l)l992f) . Using a theorem bv lKiefer and Wolfowitd l|l960|) . iFbrd et all 
l)l992j) show for a class of response curves that the D-optimal design is a two-point design. 
The proof for the complementary log-log regression is obtained as a special case of the 
general results. 



4. Sequential designs for switching measurements 

Our goal is to design a fast automated estimation procedure for switching measurements. A 
fast estimation procedure allows, for example, the tracking of changes in the parameters as 
a function of changes in environmental variables. The speed of the procedure is measured in 
terms of the total time used for measuring and estimation. Sequential designs are a natural 
choice because it is easy to control the current during the experiment. In sequential designs 
(also known as two-stage and multi-stage designs), we start with initial estimates, collect 
additional data, obtain new estimates and use them as the initial estimates for the next 
stage. 

In switching measurements, the number of stages can be high: this is a significant differ- 
ence compared to medical dose response experiments where practical reasons usually limit 
the number of stages. However, it is not necessarily advisable to choose the number of 
stages as high as possible, i.e., update the estimates after every two observations. This 
is because adding two new observations would have a very small effect on the estimates 
(especially if the number of observations is already large) but the estimation itself would 
spend time that could be used for measuring. Additionally, whenever the height of the 
current pulse is changed, the reconfiguration of the signal generator also takes a time that 
is considerably longer than the time needed for measuring the response for a single pulse. 
The optimal updating interval depends also on the computational resources and facilities, 
and our experiences suggest that it is always reasonable to measure at least 20 times using 
the same pulse height. Our recommendation is only the lower limit and based on empiri- 
cally experience on the specific measurement instruments we are using. In the experiment 
presented in Sectional we measured 50 times using the same pulse height at the beginning 
of the experiment and increased this number by 10 % at every stage so that at the stage 50 
we measured 5379 times using the same pulse height. 

Obtaining good initial estimates for the parameters is a challenging problem for switch- 
ing measurements. Naturally, the results from the previous experiments can be used to 
obtain initial estimates for the parameters but there is no guarantee that this would lead to 
good initial estimates. T he problem of poor initi al est imates in the case of logistic regres- 
sion is considered e.g. bv iKing and Wonel ((2000) and ISitterl lll992Ti . The authors propose 
sophisticated minimax approaches that have good optimality properties even if the initial 
estimates are poor. These approaches are probably applicable also to our problem but 
because of their computational complexity, we choose a simpler approach that suits espe- 
cially for switching measurements. It is usually possible to specify an interval from where 
the covariates should be taken but the problem is that this interval might be very wide 
compared to the interval where the probability of response 1 increases, for instance, from 
0.01 to 0.99. This implies that initial covariates chosen randomly from the interval may be 
very inefficient. In the worst case, poorly chosen initial covariates may lead to the situation 
where maximum likelihood estimates (MLEs) do not exist for the initial data. The existence 
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of MLEs for the model {T3J) requires that l) Albert and Andersonl Il984j) 

max(xi\yi = 0) > min(xj|j/i = 1) and 

max(xi\yi = 1) > min(xj|y, = 0). (32) 

This condition is fulfilled, for instance, if we have two points where both 0's and l's are 
measured as responses. Our approach for the initial estimation aims to finding two such 
points as quickly as possible. This heuristic procedure is can be presented as follows: 

f. Use the previous knowledge to construct such an interval [x m in,£max] that we can be 
sure that F(:r max ; a, b) — -F(x m i n ; a, b) is close to f . Constructing this kind of interval is 
usually possible even if we have very little knowledge on the parameters. 

2. Use binary search to find such a point x from the interval [x m in, x ma x] that both 0's and 
l's are measured as responses. In binary search, we measure the responses at the middle 
point of the current interval. If only 0's are measured, the middle point is taken as the 
new starting point of the interval. If only l's are measured, the middle point is taken 
as the new end point of the interval. Let x be found as the middle point of the inteval 
[xi,x u ],i.e. x = (xi + x u )/2. 

3. Measure at the point x + e, where |e| = (x u — x;)/4. The sign of e is determined according 
to the measured response: sign(e) = sign(0.5 — y), where y is the average response for x. 
If 0.5 — y = the sign of e is chosen randomly. 

4. If MLEs exist for the data measured so far, proceed to the maximum likelihood estima- 
tion. Otherwise, divide e by two and return to Step 3. 

Step 2 requires that we measure at least two times at each point. In practice, it is preferable 
to measure at least tens of times at each point because of time needed for the reconfiguration 
of the signal generator. When the binary search in Step 2 converges we have found one of 
the two points required for the maximum likelihood estimation. The other point that is 
needed is then found in the neighbourhood of point x using again binary search. Due to 
use of binary search in Steps 2 and 3, the procedure works reliably and efficiently. At every 
step of binary search, the length of the current interval is divided by two. Consequently, 
the procedure converges exponentially fast regardless of the choice of the initial interval. 

After the heuristic procedure has converged the experiment continues sequentially: The 
current values of a and b are used to determine the optimal covariates for the next stage, 
measurements are performed, and finally maximum likelihood estimation is applied to up- 
date values of a and b. In switching measurements, one is also interested in parameters that 
directly describe the location and the scale of the response curve in addition to parameters 
a and b. Let us define the middle point 9 and width of the curve A, as follows 

9=-(F- 1 (0.5)-b) (33) 
a 

\=-(F- 1 (0.9)-F- 1 (0.l)). (34) 
a 

These two parameters are frequently used in switching measurements to characterize the 
response curve. The covariance matrix of (9 A) T is given by 



(35) 
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where 



A = 




06 
db 

ax 
~d!> 




(36) 



It is often practical to define the stopping rule for the data acquisition using the variances 
of 6 and A. Special attention should be paid to the stability of the measurement environment. 
Instability, such as a change in the temperature during the experiment, causes an increase 
in the variances. This emphasizes the need for fast and reliable experimental design. 

5. Examples with real and simulated data 

We measured a sample consisting of an aluminium-aluminium oxide-aluminium JJ circuit 
in a dilution refrigerator at 20 mK temperature. The sample parameters a and b depend 
on the area of the JJ and the thickness of the oxide layer. The sample fabricated involves 
a JJ whose area is 1 times 2 micrometers and the oxide layer is ten nanometres thick. 
The sample was connected to computer controlled measurement electronics in order to 
apply the current pulses and record the resulting voltage pulses. The resistance of the 
sample at room temperature suggests that a pulse of 300 nA always causes a switching 
(response 1). This gives the upper limit for the initial estimation. The lower limit for the 
initial estimation, 200 nA, can be roughly estimated from the dimensions of the J J. The 
sequential design presented in Section 21 was used to construct the applied current pulse 
sequences in real time. Each stage consisted of measuring phase and estimation phase. 
Since the reconfiguration of our signal generator takes a relatively long time compared to 
the length of a single pulse, we used sequences of tens or hundreds of pulses to collect data at 
each stage. In the setup the pulse lengths are on the order of milliseconds. In order to keep 
the time needed for computations much below the actual measurement time, the length of 
the pulse sequence should grow exponentially. We started with a sequence of 50 pulses per 
point and let the number of pulses grow 10 percent at each stage. In the test experiments, 
the practically sufficient accuracy of the estimates was reached after 20-30 stages, implying 
that 6000-17000 pulses were needed in total. In our setup, the experiment (50 stages) took 
8 minutes 16 seconds from which 1 minute 30 seconds were used for estimation. 

Table ^ an d Fig. 0] illustrate the estimation of the unknown parameters in the exper- 
iment. In the lower panel of Fig. ^ the functional form of the response curve is verified 
empirically. The response curve is compared with individual response probabilities that are 
independently estimated using 2000 pulses per point. The results confirm that the response 
curve is given by the Gumbel distribution as suggested by the theory. 

The performance of the sequential design is also studied in a simulation example. The 
simulation aims to imitate the real data example . The response variable Y is generated 
from the model i|13|) with parameters a = 0.24 and b — —61 which are motivated by the 
actual experiment presented above. Our goal is to estimate these parameters when we ini- 
tially know only that at point £ min = 200 the response is with high probability and at 
point x max = 300 the response is 1 with high probability. The estimation starts using the 
heuristic procedure described in Section^and continues as maximum likelihood estimation. 
In the initial estimation, we measure 25 times at each point. In the sequential maximum 
likelihood estimation, the parameter estimates are updated when the total number of ob- 
servations has increased by 20 % from the previous update. The simulation is repeated 
500 times and maximum likelihood estimates d and b arc recorded for different numbers of 
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Table 1 . Estimates and their standard deviations for the mea- 
sured JJ. 



Stage 


n 




a 


std(a) 


t> 


std(6) 


1 


100 




No MLE 


No MLE 


2 


210 




No MLE 


No MLE 


3 


332 





.3800 


0.1933 


-95.5966 


48.4418 


4 


466 





.2585 


0.0392 


-65.1942 


9.8369 


5 


614 





.2301 


0.0224 


-58.1248 


5.6501 


6 


776 


0. 


,2520 


0.0177 


-63.6371 


4.4625 


7 


954 


0. 


.2474 


0.0148 


-62.5090 


3.7268 


8 


1150 


0. 


,2368 


0.0126 


-59.8039 


3.1855 


9 


1366 





2366 


0.0112 


-59.7514 


2.8222 


10 


1604 





,2426 


0.0102 


-61.2916 


2.5744 


15 


3200 





,2350 


0.0067 


-59.3476 


1.7022 


20 


5762 





,2427 


0.0050 


-61.2986 


1.2637 


25 


9890 





,2433 


0.0038 


-61.4761 


0.9572 


30 


16550 





,2424 


0.0029 


-61.2533 


0.7359 


40 


44578 


0. 


.2410 


0.0017 


-60.8778 


0.4430 


50 


117288 





,2402 


0.0011 


-60.6282 


0.2708 


Stage 


n 




6 


std(0) 


A 


std(A) 


1 


100 




No MLE 


No MLE 


2 


210 




No MLE 


No MLE 


3 


332 


250. 


,5738 


0.2519 


8.1158 


4.1282 


4 


466 


250. 


.8324 


0.3037 


11.9342 


1.8080 


5 


614 


251. 


0506 


0.2984 


13.4066 


1.3070 


6 


776 


251. 


.0632 


0.2519 


12.2392 


0.8591 


7 


954 


251. 


.1471 


0.2327 


12.4655 


0.7432 


8 


1150 


251. 


.0128 


0.2219 


13.0259 


0.6933 


9 


1366 


250. 


.9891 


0.2058 


13.0362 


0.6148 


10 


1604 


251. 


.1269 


0.1874 


12.7136 
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Fig. 4. Estimated response curve for the measured JJ. Upper panel: estimated response curve with 
data used in the estimation (pulse heights, I, and corresponding empirical probabilities, P). Lower 
panel: estimated response curve and verification points obtained as an average of responses of 
2000 independent pulses. 
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observations n = 50, n = 100, . . . ,n — 20000 . The sequential design is compared to the the- 
oretically D-optimal design where the optimal covariate values x* — (z* + 61)/0.24 = 258.25 
and X2 — (^2 + 61)/0.24 = 248.59 are assumed to be known already at the beginning of 
the experiment. As a performance measure we calculate the mean squared errors (MSE) 
between the estimates and the true parameter values. 

Table |21 summarizes the obtained MLEs and their standard deviations. The reported 
numbers are means from 500 simulations. The MSEs between the estimates and the true 
parameter values are reported in Table The MSEs are calculated only from those sim- 
ulation runs where MLEs exist. Naturally, the theoretically optimal design results smaller 
MSE values than the sequential design but it can be also seen that for large values of n 
the performance of the sequential design approaches to the performance of the theoretically 
optimal design. The results demonstrate the practical importance of good experimental 
design and the significance of reliable estimation of initial parameters. We conclude that 
the proposed sequential design performed reasonably well compared to the theoretically 
D-optimal design when n > 200. In switching measurements, the number of observations is 
always in thousands. Comparison of estimates and their standard deviations in Tableland 
Table[2]does not reveal any major differences between the experimental data and simulated 
data. This supports the practical usability of the sequential design. 

6. Conclusion 

We have considered optimal designs for switching measurements of superconducting Joseph- 
son junction circuits. The D-optimality criterion suggests that the measurements at points 
where ax\ + b = 0.979633 and ax2 + b = —1.33774 will yield a maximal information about 
the estimated parameters assuming that the probability distribution of the measurement 
values are given by Eq. As a practical method for switching measurements, we propose 
a sequential design that combines heuristic search for initial estimates and maximum like- 
lihood estimation. The fast data acquisition allows tracking of fast environmental changes. 
The presented experimental design is now in daily use in Low Temperature Laboratory at 
Helsinki University of Technology. 

The current work opens several directions for the future research. Besides D-optimality 
we may also consider other optimality criteria. These criteria likely lead to designs different 
from the D-optimal design but because of the complicated functional forms it may be 
difficult to derive analytical results. 

The problem of choosing the number of measurements per stage was approached em- 
pirically in this paper. A more systematic approach would determine the optimal number 
of measurements per stage as a function of the observed data. It seems reasonable to as- 
sume that each measurement takes a certain time and the change of the measurement point 
takes additional fixed time. The problem is to choose the number of measurements per 
stage in such a way that the expected total time needed to achieve the required accuracy 
is minimized. This is a non-trivial problem and probably only approximate solutions are 
reachable. 

The time used for the experiment could be further reduced if the estimation and the 
measuring were done simultaneously. This requires good coordination of data processing and 
knowledge on the time needed for each operation. In sequential designs, the simultaneous 
estimation and measuring implies that the current estimates are not based on all data 
measured so far but all data processed so far. 
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Table 2. Means of estimates and their standard deviations from 500 simulation runs. 
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Table 3. Mean squared errors (MSE) for parameter estimates 
a and b as functions of number of observations n. The proposed 
sequential design and the theoretically D-optimal design are com- 
pared. Reported MSEs are means from 500 simulation runs. 
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0.399793 



Although the presented complementary log-log regression model describes well some JJ 
circuits involving several junctions, we probably need more flexible models for complicated 
structures of JJ. The modelling can be approached from two different directions: we may 
derive the exact physical model and try to find the optimal design for it, or we may start 
with a flexible parametric model for which the optimal design is known and estimate the 
model parameters from the measured data. 

Statistical methods can be also applied to the modelling of JJ parameters as a function 
of environmental variables, such as temperature. This direction of research leads to adaptive 
sequential designs where the challenge is not only to estimate the parameters but also to 
track their changes in optimal way. 
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